Nuclear quantum effects in water 
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In this work, a path integral Car-Parrinello molecular dynamics simulation of liquid water is 
performed. It is found that the inclusion of nuclear quantum effects systematically improves the 
agreement of first principles simulations of liquid water with experiment. In addition, the proton 
momentum distribution is computed utilizing a recently developed open path integral molecular 
dynamics methodology. It is shown that these results are in good agreement with neutron Compton 
scattering data for liquid water and ice. 

PACS numbers: 61.20.Ja, 71.15.Pd 
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Due to the fundamental importance of water in the 
physical and biological sciences, understanding its mi- 
croscopic structure is an issue of longstanding interest. 
Elucidating the local environment of the protons is par- 
ticularly intriguing due to their crucial role in hydrogen 
bonding. Nuclear quantum effects have a significant im- 
pact on the behavior of water. This is indicated by the 
large isotope effects observed in numerous water prop- 
erties when protons are substituted with deuterium (D) 
or tritium (T) nuclei. For example, the melting point of 
heavy water (D 2 0) is 3.82 K higher than that of regular 
(H 2 0) water, and the effect is even more pronounced in 
tritiated water (T 2 0) fl, providing evidence that quan- 
tum effects destabilize the hydrogen bond network. 

Recently, the equilibrium state of the protons in water 
and ice has beenprobed by neutron Compton scatter- 
ing experiments [2| . This technique measures the proton 
momentum distribution [3J], thereby providing comple- 
mentary information to what is garnered from diffraction 
techniques that measure the spatial correlations among 
the nuclear positions 0, H| . Due to the non-commuting 
character of position and momentum operators in quan- 
tum mechanics, the proton momentum distribution is 
sensitive to the local environment. In particular, the dif- 
ferences in the momentum distribution of the solid and 
liquid water phases reflect the breaking and distortion of 
hydrogen bonds that occurs upon melting. In systems 
such as confined water @,J3 and the quantum ferroelec- 
tric potassium phosphate [8j , the momentum distribution 
provides signatures of tunneling and derealization. 

Molecular simulations with quantum nuclei are made 
feasible by the Feynman path integral representation of 
the equilibrium density matrix at finite temperature. 
This approach has been used in conjunction with em- 
pirical force fields in studies [j| showing that quantum 
fluctuations soften the structure of liquid water. The ef- 



fect is illustrated by a broadening of the radial distribu- 
tion functions (RDF) compared to those generated from 
classical nuclei. Interestingly, these works indicated that 
quantum nuclei affect the structure in a similar way to a 
temperature increase in a classical simulation. Recently, 
empirical force fields have been employed within "open" 
path integral molecular dynamics methodologies to com- 
pute the proton momentum distribution in ice and wa- 
ter 13, 11 1 . The calculated distribution, while in agree- 
ment with experiment in many respects, did not repro- 
duce the shorter tail that is observed in ice, signaling a 
lack of transferability of the empirical potentials. The 
faster decaying ice distribution reflects a red-shift of the 
OH stretch frequency that is a consequence of the recov- 
ery of an intact hydrogen bond network upon freezing. 

To investigate whether this effect can be reproduced 
in ab-initio simulations, we perform an "open" path inte- 
gral Car-Parrinello molecular dynamics (PI-CPMD) [l2| 
study of water in the liquid and solid phases. In this 
approach the nuclear potential energy surface is de- 
rived on the fly from the instantaneous ground state of 
the electrons within Density Functional Theory (DFT). 
Our study is also motivated by a previous, pioneer- 
ing PI-CPMD simulation of liquid water [l3j. This 
study reached the counterintuitive conclusion that nu- 
clear quantum effects harden the structure of the liquid 
in comparison to classical CPMD simulations at the same 
temperature. Numerous studies have shown that such 
simulations generate an overstructured liquid 
Consequently, nuclear quantum effects would increase the 
discrepancy between experiment and simulation. If cor- 
rect, this result would have severe implications for the 
accuracy of current DFT approximations of water. 

In this work we use a combination of closed and open 
Feynman paths to compute the pair correlation functions 
and the momentum distribution. We find that the liquid 
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is significantly less structured than in computations uti- 
lizing an identical electronic structure description with 
classical nuclei, in qualitative agreement with experi- 
mental isotope effects and previous force field studies. 
The computed proton momentum distributions are also 
in good agreement with experiment and, unlike in em- 
pirical force field based studies, the difference between 
the liquid and the solid observed in experiment is repro- 
duced. Small remaining deviations from experiment sug- 
gest some degree of over-binding in the hydrogen bond 
network that is likely engendered by the adopted approx- 
imate DFT description of the electronic structure. 

In the primitive discretization of the path integral for- 
malism, the problem of describing the quantum nuclei is 
mapped onto a set of classical replicas coupled via har- 
monic interactions. In order to compute properties such 
as the momentum distribution that are not diagonal in 
the position representation, we must use an open path to 
represent the nuclei for which the momentum distribu- 
tion is computed [lTf. All other nuclei should be repre- 
sented by closed paths. The proton momentum distribu- 
tion may then be computed from the Fourier transform 
of the open path end-to-end distribution. In the present 
study 32 replicas are employed 1(J II |- We utilize the 
algorithm presented in Ref. [O to perform the nuclear 
dynamics. In this algorithm, one hydrogen per water 
molecule is represented by an open path, an approxi- 
mation that facilitates efficient sampling with insignif- 
icant impact upon the momentum distribution. This 
scheme was implemented within the CPMD [18[ pack- 
age, which has been optimized for the IBM Blue Gene/L 
platform on which these simulations were carried 
out. The staging transformation [2(| and massive Nose- 
Hoover thermostat chains [2ll ] are employed in order to 
ensure the sustained diffusion of the trajectory, indicated 
by the mean square displacement of the oxygen centroid. 
The fictitious sampling masses are set to be 4 times larger 
than the corresponding physical masses of the nuclei. 

The atomic forces are computed from first principles 
via the CPMD methodology [2^]. Electron exchange 
and correlation effects are treated with the BLYP func- 
tional 2^, 2^ 1 . T roullier-Martins norm-conserving pseu- 
dopotentials 25| are employed and the Kohn-Sham or- 
bitals are expanded in a plane wave basis set with a 75 
Ry cutoff. A fictitous electron mass of 340 atomic units 
and a time step of 3.0 atomic units is chosen [26]. 

The liquid water system contains 64 molecules simu- 
lated at 300K and is placed in a periodic cubic box of 
length 12.459A, yielding a density within 0.6% of the ex- 
perimental value at the simulation temperature. After 
an equilibration period of 6 ps, a production run of 12.6 
ps is generated. An electron thermostat is employed in 
order to maintain the fictitious electronic kinetic energy 
at an average value of 32K per nuclear degree of free- 
dom. In order to assess the effect of quantum nuclei, a 
CPMD simulation of liquid water with classical nuclei is 



also carried out utilizing otherwise identical parameters 
at temperatures of 300K and 330K, with average fictitous 
electronic kinetic energies of 17K and 18K, respectively. 
A 21 ps production trajectory is generated after an 8 ps 
equilibration period. The proton-disordered hexagonal 
ice system contains 96 molecules simulated at 269K and 
is modeled at experimental density in a periodic box with 
one side of length 13.556A and the sides in a 1 : 1.15 : 1.09 
ratio. The electronic kinetic energy is held to a value of 
27K. After 1 ps of equilibration, a 3.8 ps production run 
is generated. 

The 00 and OH RDFs for liquid water are shown in 
Fig. [TJ They are plotted against RDFs garnered from 
neutron scattering experiments The inclusion of nu- 
clear quantum effects leads to significantly less struc- 
tured RDFs than the corresponding CPMD simulation 
with classical nuclei. This qualitative observation is in 
agreement with previous studies of water that employ 
empirical potentials [lj|. The covalent peak of the OH 
RDF is slightly shifted, but in otherwise good agreement 
with the experimental neutron scattering results, though 
the second peak, which corresponds to hydrogen bonding 
interactions, is somewhat sharper, although not nearly as 
sharp, as the standard CPMD result. The first peak of 
the OO RDF is a useful marker of the relative structuring 
of water. The resultant peak from the path integral sim- 
ulation is 2.84, which is much closer to thepeak height 
garnered from both neutron 0] and x-ray Q scattering 
experiments than the peak height of 3.23 that we obtain 
from a standard CPMD simulation. The latter result is 
in excellent agreement with a previous study that em- 
ployed a similar methodology and parameters [l6j |. The 
path integral 00 RDF is overall very similar to that of 
a standard simulation run at 330K. However, the quan- 
tum derealization of the protons is stronger than the 
classical temperature effect at 330K, as indicated by the 
OH RDFs in Fig. [1] From these results, it appears that 
the overstructuring present in standard CPMD simula- 
tions at room temperature is in part mitigated by the 
inclusion of nuclear quantum effects. Yet, despite the im- 
provement, there remains some degree of overstructuring 
in the path integral result. In particular, the first min- 
imum of the 00 RDF is deeper than the experimental 
result. 

One may also evaluate the degree of structuring in 
the liquid from the average fraction of broken hydrogen 
bonds. A hydrogen bond is defined in geometric terms by 
oxygen-to-oxygen and oxygen-to-hydrogen distance cut- 
offs that are equal to the minima of the hydrogen-bonding 
peaks of the RDFs (Fig. [T]), and a hydrogen bond angle 
greater than 140 degrees. In the path integral represen- 
tation, the distances and angles are measured from the 
centroid of each representative path. Absence of broken 
hydrogen bonds yields a tetrahedral coordination for each 
molecule as in ice. The fraction is approximately 10% for 
liquid water [5j, and both the path integral (11%) and 
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FIG. 1: (color online) The OO (top panel) and OH (bottom 
panel) RDFs in liquid water from an open PI-CPMD simula- 
tion at 300K (solid line), standard CPMD simulations at 300K 
(dashed line) and 330K (dotted line), and neutron scattering 
experiments at 298K 4] (circles) are reported. 

standard simulation (7%) yield averages near this value. 
However, the larger number of broken hydrogen bonds in 
the PI-CPMD result indicate increased fluidity. 

The distribution of dipole moments in the standard 
and path integral simulation at 300K is depicted in Fig. 
[2J The dipole moment of each individual molecule was 
compiled over selected configurations via the sum over 
ions and the centers of maximally localized Wannier func- 
tions [13]. Within the available statistics there is no ap- 
preciable change in the average dipole moment of the 
classical and quantal distributions. The only noticeable 
difference is a broadening of the path integral distribution 
due to the broadening of the OH covalent bond distribu- 
tion evident in Fig. Q] As a consequence, the root-mean 
square dipole moment is larger in the path integral than 
in the standard simulation. Since the dielectric constant 
is proportional to the fluctuation of the molecular dipole 
moment [28j , this result is consistent with the experimen- 
tal finding that the dielectric constant of light water is 
slightly larger than that of heavy water [lj . 

It is important to assess how the equilibrium properties 
of position are affected by our "open" path approxima- 
tion. The OH RDF is, when calculated with open path 
hydrogens, broader by approximately 4% than when cal- 
culated with closed path hydrogens. However, if the four 
end replicas on each side of an open chain are excluded 
from the average, the corrected open OH distribution is 
similar to the one corresponding to closed paths, which 
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FIG. 2: The dipole moment distribution of the full PI simu- 
lation (solid line), corrected PI simulation (dotted line) and 
standard (dashed line) simulation of water at 300K. 

is reported in Fig. Q] Such effects are less pronounced 
for the OO radial distribution function. Open and closed 
hydrogen chains are used to calculate the dipole moment 
distribution in Fig. [2| where a corrected plot is also re- 
ported. One notes that the corresponding distribution is 
essentially the same as the uncorrected one. We also find 
that the open path approximation has a negligible effect 
on the computed fraction of broken hydrogen bonds. 

The proton momentum distribution is computed for 
liquid water and ice, and is compared to neutron Comp- 
ton scattering data [2| and the Boltzmann distribution in 
Fig. [3l There is a large distinction between the classical 
result, which only depends upon mass and temperature, 
and the actual momentum distribution which, owing to 
quantum effects, is sensitive to the potential energy sur- 
face. In both phases, the proton momentum distribution 
of the simulation is broader than the experimental result, 
although the curves are generally in good agreement with 
each other. Consistent with the uncertainty relation be- 
tween position and momentum, the broader computed 
momentum distribution corresponds to the more struc- 
tured (i.e. more localized) OH RDF of the simulation in 
comparison to experiment. This is depicted for the liquid 
in Fig. [I] and is also noticed upon inspection of the ice 
OH RDF (not shown). The tail of the computed distribu- 
tion in ice is shorter than in the liquid, which can be seen 
from the insets of Fig. [3l and is in good qualitative agree- 
ment with the experimental distributions. This effect has 
not been reproduced in open path i nteg ral simulations 
that employ empirical force fields 0, [II]] . Therefore the 
first principles potential energy surface, unlike less trans- 
ferable interaction models, is sensitive to the red-shift in 
the OH stretch of a water molecule when it is placed in 
a stronger hydrogen bonding environment. 

Our position and momentum results are mutually con- 
sistent and agree with both physical intuition and the 
available experimental data. We find that nuclear quan- 
tum effects considerably soften the structure of the liquid 
and consequentially correct in large part the overstruc- 
turing present in standard first principles water. Even 
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FIG. 3: (color online) The radial proton momentum distribu- 
tions is reported in the liquid (solid line) and solid (dashed 
line) phases and plotted against the experimental liquid (cir- 
cles) and solid (triangles) curves 0] , as well as the Boltzmann 
distribution (dashed line) at 300K. The insets (a) and (b) de- 
pict the tail of the liquid (solid line) and solid (dashed line) 
distributions, in the simulation and experiment, respectively. 

with the inclusion of such effects, however, there is still 
some degree of overstructuring present in the simulation, 
indicating a residual error in the description of the po- 
tential energy surface. This is likely due to the adopted 
treatment of exchange and correlation, as studies have 
shown that hybrid functionals improve the description of 
hydrogen bonding [29(. It has also been suggested that 
the use of plane wave basis sets at typical cutoffs con- 
tributes to the overstructuring of liquid water [30| . 

In this work, we have extended the PI-CPMD method- 
ology to allow computations of the proton momentum 
distribution and have reported the first application of 
this scheme to liquid and solid water. Our results are in 
good agreement with neutron Compton scattering data. 
Given the similarity of the local environment in water and 
ice, the improvement provided by this approach over em- 
pirical potentials, albeit qualitatively important, is quan- 
titatively modest. This approach should be particularly 
useful in treating proton tunnelling, which involves bond 
breaking and forming events that are not easily captured 
by empirical potentials. Such events are likely to play a 
crucial role in proton wires that are present in biological 
settings [3l|. Contributing to the understanding of such 
systems is a future goal of both this methodology and 
the related experimental techniques. 
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of their computational resources. 
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